
import numpy as np
import matplotlib.pyplot as plot
import math
import random
def get_bud(a,b,n):
     bud = (b-a) / (n-1)
     return bud
#---------------相关参数-------------------------
a = -10
b = 10
c = 2
d = 3
e = 5
f = 4
n = 100
m = 10


#构建插值函数
bud = get_bud(a,b,n)
scalex = np.arange(a,b+bud,bud)
bud2 = get_bud(a,b,m)
ylist =[]
for i in scalex:
    ylist.append(c*math.sin(d*i) + e * math.cos(f*i))


#构建测试点
#在区间内取m个均匀的点
pf =a;
x_test = [ i for i in range(0,m)]
for i in range(0,m):
    x_test[i] = pf
    pf += bud2


y_test = []
for i in x_test:
    y_test.append(c*math.sin(d*i) + e * math.cos(f*i))



#最小二乘法函数拟合
#x[]为离散点集合的x值，y[]为离散点集合的y值
#s为最高项次数
def get_p(x,fx,s):
    row = len(x)
    low = s + 2
    p = [[0 for i in range(low)]for j in range(row)]
    #p[row][low] 
    #p[row][0]存放x的值
    #p[row][1]表示多项式0次幂时的系数
    for i in range(0,row):
        p[i][0] = x[i]
        p[i][1] = 1
    #对每一个x进行计算ak和bk并计算Pk
    for i in range(2,low):
        a1=0
        a2=0
        b1=0
        b2=0
        a = [i for i in range(low)]
        b = [i for i in range(low)]
        for j  in range(0,row):
            for m in range(0,row):
               a1 += 1 * p[m][0] * p[m][i-1] * p[m][i-1]
               a2 += 1 * p[m][i-1] * p[m][i-1]
               if(i >=3):
                 b1 += p[m][i-1] * p[m][i-1]
                 b2 += p[m][i-2] * p[m][i-2]
               else:
                b[i-2] = 0   
            if(i>=3):
              b[i-2] = b1 / b2  
            a[i-1] =  a1 / a2
            p[j][i] = (x[j] - a[i-1]) * p[j][i-1]  - b[i-2] * p[j][i-2] 
    A = [i for i in range(low)]
    for i in range(1,s+2):
        A1 = 0
        A2 = 0
        for j in range(0,row):
            A1+= 1 * (fx[j] ) * p[j][i]
            A2+= 1* p[j][i] * p[j][i]
        A[i] = A1 / A2
    y = [i for i in range(row)]
    for i in range(row):
        y[i] = 0
    for j in range(0,row):
       for i in range(1,s+2):
           y[j] += A[i] * p[j][i]
    return y

row = len(scalex)
low = 4
y1 = get_p(scalex,ylist,40)
#计算误差
y_list = get_p(x_test,y_test,40)
sum = 0
for i in range(0,len(x_test)):
    sum+=(y_test[i] - y_list[i])
print(sum / len(x_test))
plot.plot(scalex,y1,color = 'blue',label = 'oridinary')
plot.plot(scalex,ylist,'ro')
plot.legend()
plot.show()

